Python with PostgreSQL & PostGIS¶

Note: Please always run the complete Jupyter Notebook from the beginning, as object names such as 'sql' and 'gdf' are reused in the code cells.

Libraries and Settings¶

The 3 types of swiss coordinate systems are

We adapt the code in this part: ST_transform(p.way, 4326) as geom

We have either 21781, 20256 or 4326

This where I can find the information regarding shops code

https://wiki.openstreetmap.org/wiki/Key:shop

TO get to coordinates https://tools.retorte.ch/map/

In [1]:
# Libraries
import os
import folium
import pandas as pd
import geopandas as gpd
from sqlalchemy import create_engine, text

# Ignore warnings
import warnings
warnings.filterwarnings("ignore")

print(os.getcwd())

/workspaces/python_postgresql_postgis

Create database connection¶

In [2]:
# Set up database connection
user = "pgadmin"
password = "geheim"
host = "localhost"
port = "5432"
database = "osm_switzerland"

# Create Connection URL
db_connection_url = "postgresql://" + user + ":" + password +\
                    "@" + host + ":" + port + "/" + database

# Create SQLAlchemy Engine
engine = create_engine(db_connection_url)

# Test database connection
with engine.connect() as connection:
    result = connection.execute(text('SELECT current_database()'))
    print(result.fetchone())

# Dispose the engine
engine.dispose()

('osm_switzerland',)

List tables in database¶

In [3]:
# Create SQLAlchemy Engine
engine = create_engine(db_connection_url)

# Open a connection
with engine.connect() as connection:

    # Execute the query
    result = connection.execute(text("""SELECT table_name
                                        FROM information_schema.tables
                                        WHERE table_schema = 'public';"""))
    
    # Fetch and print the results
    for row in result:
        print(row[0])

# Dispose the engine
engine.dispose()

geography_columns
geometry_columns
spatial_ref_sys
planet_osm_roads
planet_osm_line
planet_osm_polygon
planet_osm_point

Show columns and data types of selected table¶

In [4]:
# Create SQLAlchemy Engine
engine = create_engine(db_connection_url)

# Specify your table name
table_name = 'planet_osm_polygon'

# Query to get column information
query = f"""SELECT column_name, data_type 
        FROM information_schema.columns 
        WHERE table_name = '{table_name}'"""

# Execute the query and read the result into a DataFrame
df = pd.read_sql(query, engine)

# Dispose the engine
engine.dispose()

# Print the DataFrame
df
Out[4]:
column_name data_type
0 osm_id bigint
1 z_order integer
2 way_area real
3 way USER-DEFINED
4 addr:housenumber text
... ... ...
68 wood text
69 tracktype text
70 access text
71 addr:housename text
72 addr:street text

73 rows × 2 columns

Query: Select buildings for which full address is available in defined zip code areas¶

In [5]:
# Create SQLAlchemy Engine
engine = create_engine(db_connection_url)

# Define SQL query 
sql = """SELECT
                p.osm_id,
                p."addr:street",
                p."addr:housenumber",
                p."addr:city",
                p."addr:postcode",
                p.building,
                st_transform(p.way, 4326) AS geom
        FROM
                public.planet_osm_polygon AS p
        WHERE
                p."addr:street" IS NOT NULL
                AND p."addr:housenumber" IS NOT NULL
                AND p."addr:city" IS NOT NULL
                AND p."addr:postcode" IN ('8001', '8002')"""

# Query the database and store the result in a GeoDataFrame
gdf = gpd.GeoDataFrame.from_postgis(sql, engine)

# Dispose the engine
engine.dispose()

# Print the GeoDataFrame
gdf
Out[5]:
osm_id addr:street addr:housenumber addr:city addr:postcode building geom
0 39517505 Schanzengasse 14 Zürich 8001 yes POLYGON ((8.54989 47.36681, 8.54996 47.36676, ...
1 288353536 Fortunagasse 34 Zürich 8001 yes POLYGON ((8.53967 47.37337, 8.5397 47.37334, 8...
2 288353533 Fortunagasse 30 Zürich 8001 yes POLYGON ((8.53967 47.37346, 8.53969 47.37344, ...
3 288353542 Rennweg 42 Zürich 8001 yes POLYGON ((8.53931 47.37351, 8.53933 47.37347, ...
4 288353540 Rennweg 38 Zürich 8001 yes POLYGON ((8.53933 47.37347, 8.53939 47.37338, ...
... ... ... ... ... ... ... ...
2350 108633058 Brandschenkestrasse 152 Zürich 8002 industrial POLYGON ((8.52413 47.36489, 8.52415 47.36477, ...
2351 108633096 Brandschenkestrasse 152c Zürich 8002 office POLYGON ((8.52388 47.36409, 8.52391 47.3639, 8...
2352 108633055 Brandschenkestrasse 152b Zürich 8002 industrial POLYGON ((8.5239 47.36435, 8.5239 47.36431, 8....
2353 108626600 Brandschenkestrasse 110a Zürich 8002 office POLYGON ((8.52458 47.3655, 8.52462 47.36534, 8...
2354 34572240 Brandschenkestrasse 110 Zürich 8002 commercial POLYGON ((8.52469 47.36579, 8.52476 47.36552, ...

2355 rows × 7 columns

Show selected features on map¶

Note the popup field in the map, which has been added to provide additional information about buildings.

Example of alternative background maps (maptiles) are:

  • EsriWorldImagery
  • EsriWorldTopoMap
  • EsriWorldGrayCanvas
  • CartoDBDarkMatter
  • CartoDBPositron
In [6]:
# Ensure the GeoDataFrame has the correct projection (EPSG:4326)
if gdf.crs is None:
    gdf.set_crs(epsg=4326, inplace=True)
else:
    pass

# Calculate the mean longitude and latitude for the map center
centroids = gdf.geometry.centroid
lon = centroids.x.mean()
lat = centroids.y.mean()

# Initialize the map
m = folium.Map(location=[lat, lon], 
               zoom_start=15,
               tiles='EsriWorldImagery')

# Map settings
folium.GeoJson(
    gdf,
    name='geojson',
    weight=0.5,
    fill_color='greenyellow',
    fillOpacity=0.8,
    popup=folium.GeoJsonPopup(fields=['addr:street',
                                      'addr:housenumber',
                                      'addr:city',
                                      'addr:postcode',
                                      'building'])
).add_to(m)

folium.LayerControl().add_to(m)

# Plot map
m
Out[6]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Doing the same exercise for my city: Treytorrens¶

In [7]:
# Create SQLAlchemy Engine
engine = create_engine(db_connection_url)

# Define SQL query 
sql = """SELECT
                p.osm_id,
                p."addr:street",
                p."addr:housenumber",
                p."addr:city",
                p."addr:postcode",
                p.building,
                st_transform(p.way, 4326) AS geom
        FROM
                public.planet_osm_polygon AS p
        WHERE
                p."addr:street" IS NOT NULL
                AND p."addr:housenumber" IS NOT NULL
                AND p."addr:city" IS NOT NULL
                AND p."addr:postcode" IN ('1530')"""

# Query the database and store the result in a GeoDataFrame
gdf = gpd.GeoDataFrame.from_postgis(sql, engine)

# Dispose the engine
engine.dispose()

# Print the GeoDataFrame
gdf
Out[7]:
osm_id addr:street addr:housenumber addr:city addr:postcode building geom
0 690348871 Aéropôle 201 Payerne 1530 industrial POLYGON ((6.91756 46.83402, 6.91788 46.83383, ...
1 138739158 Grand'Rue 16 Payerne 1530 yes POLYGON ((6.93893 46.82179, 6.93904 46.82161, ...
2 138739168 Grand'Rue 20 Payerne 1530 yes POLYGON ((6.93863 46.82193, 6.93868 46.82187, ...
3 138739172 Grand'Rue 22 Payerne 1530 yes POLYGON ((6.93866 46.82174, 6.93872 46.82158, ...
4 138739160 Rue des Blanchisseuses 18 Payerne 1530 apartments POLYGON ((6.93818 46.82203, 6.93828 46.8219, 6...
... ... ... ... ... ... ... ...
1045 357477601 Rue des Jumelles 22 Payerne 1530 yes POLYGON ((6.9367 46.81495, 6.9368 46.81489, 6....
1046 357477604 Rue des Jumelles 24 Payerne 1530 yes POLYGON ((6.93651 46.81481, 6.9366 46.81474, 6...
1047 357477606 Rue des Jumelles 25 Payerne 1530 yes POLYGON ((6.93677 46.81465, 6.93685 46.8146, 6...
1048 357477610 Rue des Jumelles 29 Payerne 1530 yes POLYGON ((6.93669 46.81428, 6.9368 46.81422, 6...
1049 357477608 Rue des Jumelles 27 Payerne 1530 yes POLYGON ((6.93634 46.81433, 6.93643 46.81428, ...

1050 rows × 7 columns

In [8]:
# Ensure the GeoDataFrame has the correct projection (EPSG:4326)
if gdf.crs is None:
    gdf.set_crs(epsg=4326, inplace=True)
else:
    pass

# Calculate the mean longitude and latitude for the map center
centroids = gdf.geometry.centroid
lon = centroids.x.mean()
lat = centroids.y.mean()

# Initialize the map
m = folium.Map(location=[lat, lon], 
               zoom_start=17,
               tiles='EsriWorldTopoMap')

# Map settings
folium.GeoJson(
    gdf,
    name='geojson',
    weight=0.5,
    fill_color='red',
    fillOpacity=0.9,
    popup=folium.GeoJsonPopup(fields=['addr:street',
                                      'addr:housenumber',
                                      'addr:city',
                                      'addr:postcode',
                                      'building'
                                      ])
).add_to(m)

folium.LayerControl().add_to(m)

# Plot map
m
Out[8]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Query: Select coffee stores in Switzerland¶

In [24]:
# Create SQLAlchemy Engine
engine = create_engine(db_connection_url)  

# Define SQL query
sql = """SELECT
            h.osm_id,
            h.shop,
            h.name,
            ST_Transform(h.way, 4326) AS geom
        FROM planet_osm_point h
        WHERE h.shop = 'coffee';"""

# Query the database and store the result in a GeoDataFrame
gdf = gpd.GeoDataFrame.from_postgis(sql, engine)

# Dispose the engine
engine.dispose()

# Print the GeoDataFrame
gdf
Out[24]:
osm_id shop name geom
0 10275946038 coffee Quinta Coira POINT (9.53125 46.8485)
1 11920528580 coffee Maison du Cafe POINT (8.98798 45.87081)
2 9780147859 coffee Bel-Air Coffee POINT (6.62851 46.52281)
3 11930216622 coffee Charbon Café POINT (6.6277 46.51435)
4 7240932670 coffee Nespresso POINT (6.69463 46.50053)
... ... ... ... ...
113 11993710081 coffee Berger Kaffee POINT (7.62597 46.87853)
114 1438791023 coffee Nespresso POINT (7.63094 46.75724)
115 1517037969 coffee Nespresso POINT (7.15165 46.80227)
116 5019777080 coffee KAFOJ POINT (7.2462 47.14107)
117 5611649752 coffee Pelluch GmbH Kaffeemaschienen POINT (7.54718 47.54892)

118 rows × 4 columns

Show selected features on map¶

In [25]:
# Ensure the GeoDataFrame has the correct projection (EPSG:4326)
if gdf.crs is None:
    gdf.set_crs(epsg=4326, inplace=True)
else:
    pass

# Calculate the mean longitude and latitude for the map center
centroids = gdf.geometry.centroid
lon = centroids.x.mean()
lat = centroids.y.mean()

# Initialize the map
m = folium.Map(location=[lat, lon], 
               zoom_start=9, 
               tiles='EsriWorldTopoMap')

# Map settings
folium.GeoJson(
    gdf,
    name='map',
    popup=folium.GeoJsonPopup(fields=['name', 'shop'])
).add_to(m)

folium.LayerControl().add_to(m)

# Plot map
m
Out[25]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Same but for dairy products¶

In [11]:
# Create SQLAlchemy Engine
engine = create_engine(db_connection_url)  

# Define SQL query
sql = """SELECT
            h.osm_id,
            h.shop,
            h.name,
            ST_Transform(h.way, 4326) AS geom
        FROM planet_osm_point h
        WHERE h.shop = 'dairy';"""

# Query the database and store the result in a GeoDataFrame
gdf = gpd.GeoDataFrame.from_postgis(sql, engine)

# Dispose the engine
engine.dispose()

# Print the GeoDataFrame
gdf
Out[11]:
osm_id shop name geom
0 3086652026 dairy Käserei Sax POINT (9.45643 47.2294)
1 6331614596 dairy Nossa Caschareia POINT (9.59642 46.59588)
2 7796298536 dairy Laiterie POINT (6.48589 46.86461)
3 12090655227 dairy Laiterie POINT (6.89998 46.52691)
4 9628876022 dairy Laiterie Modèle POINT (7.01441 46.25005)
... ... ... ... ...
76 980819667 dairy Käserei Büel POINT (7.28702 46.74436)
77 4307193372 dairy Laiterie - Fromagerie Pasquier POINT (7.06689 46.76837)
78 1512803412 dairy None POINT (7.1554 47.36553)
79 10298620447 dairy Dorfkäserei Thundorf POINT (8.96682 47.54753)
80 2990133771 dairy Chäs-Hütte POINT (9.01833 47.5641)

81 rows × 4 columns

In [12]:
# Ensure the GeoDataFrame has the correct projection (EPSG:4326)
if gdf.crs is None:
    gdf.set_crs(epsg=4326, inplace=True)
else:
    pass

# Calculate the mean longitude and latitude for the map center
centroids = gdf.geometry.centroid
lon = centroids.x.mean()
lat = centroids.y.mean()

# Initialize the map
m = folium.Map(location=[lat, lon], 
               zoom_start=9, 
               tiles='CartoDBDarkMatter')

# Map settings
folium.GeoJson(
    gdf,
    name='map',
    popup=folium.GeoJsonPopup(fields=['name', 'shop'])
).add_to(m)

folium.LayerControl().add_to(m)

# Plot map
m
Out[12]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Query: Select all supermarkets in a distance of 1000m around the central station in the city of Winterthur.¶

Note:

For each supermarket, the distance to the central station in meters is calculated and stored as new column 'distance_meters'.

In addition, a popup field was added to the map, allowing users to view detailed information about each selected feature when they click on it.

The WGS84 (World Geodetic System 1984) coordinates in ST_MakePoint(LON, LAT) were derived from: https://tools.retorte.ch/map.

In [13]:
# Create SQLAlchemy Engine
engine = create_engine(db_connection_url)  

# Define SQL query
sql = """SELECT
            p.osm_id,
            p.shop,
            p.name,
            ST_Distance(
                ST_Transform(p.way, 4326)::geography,
                -- Central station coordinates
                ST_SetSRID(ST_MakePoint(8.72397, 47.50031), 4326)::geography
            ) AS distance_meters,
            ST_TRANSFORM(p.way, 4326) AS geom
        FROM
            planet_osm_point AS p
        WHERE
            p.shop = 'supermarket'
            AND ST_DWithin(
                ST_Transform(p.way, 4326)::geography,
                -- Central station coordinates
                ST_SetSRID(ST_MakePoint(8.72397, 47.50031), 4326)::geography,
                1000
            )
        ORDER BY distance_meters;"""

# Query the database and store the result in a GeoDataFrame
gdf = gpd.GeoDataFrame.from_postgis(sql, engine)

# Dispose the engine
engine.dispose()

# Print the GeoDataFrame
gdf
Out[13]:
osm_id shop name distance_meters geom
0 706203439 supermarket Coop 159.883419 POINT (8.72594 47.50085)
1 4109460421 supermarket Asia Shop 162.391281 POINT (8.72208 47.50101)
2 3831772784 supermarket Migros 247.578208 POINT (8.72115 47.49916)
3 7380954145 supermarket Alnatura 256.838011 POINT (8.72074 47.49958)
4 4095400190 supermarket ALDI 274.275393 POINT (8.72476 47.4979)
5 4125136758 supermarket Tandoor Indischer Supermarkt 290.212664 POINT (8.72017 47.50073)
6 4095400136 supermarket Denner 316.354037 POINT (8.72036 47.49886)
7 709022324 supermarket Claro Weltladen 441.129317 POINT (8.72912 47.49842)
8 4058248551 supermarket Migros 600.117307 POINT (8.73193 47.50012)
9 3441033104 supermarket L'Ultimo Bacio 680.202961 POINT (8.73299 47.49999)

Show selected features on map¶

In [14]:
# Ensure the GeoDataFrame has the correct projection (EPSG:4326)
if gdf.crs is None:
    gdf.set_crs(epsg=4326, inplace=True)
else:
    pass

# Calculate the mean longitude and latitude for the map center
centroids = gdf.geometry.centroid
lon = centroids.x.mean()
lat = centroids.y.mean()

# Initialize the map
m = folium.Map(location=[lat, lon], 
               zoom_start=16, 
               tiles='ESRIWorldImagery')

# Map settings
folium.GeoJson(
    gdf,
    name='map',
    popup=folium.GeoJsonPopup(fields=['name', 'distance_meters'])
).add_to(m)

folium.LayerControl().add_to(m)

# Plot map
m
Out[14]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Same but for coffee aroudn my home¶

In [15]:
# Create SQLAlchemy Engine
engine = create_engine(db_connection_url)  

# Define SQL query
sql = """SELECT
            p.osm_id,
            p.shop,
            p.name,
            ST_Distance(
                ST_Transform(p.way, 4326)::geography,
                -- my coordinates
                ST_SetSRID(ST_MakePoint(8.522117, 47.370249), 4326)::geography
            ) AS distance_meters,
            ST_TRANSFORM(p.way, 4326) AS geom
        FROM
            planet_osm_point AS p
        WHERE
            p.shop = 'coffee'
            AND ST_DWithin(
                ST_Transform(p.way, 4326)::geography,
                -- my coordinates
                ST_SetSRID(ST_MakePoint(8.522117, 47.370249), 4326)::geography,
                2000
            )
        ORDER BY distance_meters;"""

# Query the database and store the result in a GeoDataFrame
gdf = gpd.GeoDataFrame.from_postgis(sql, engine)

# Dispose the engine
engine.dispose()

# Print the GeoDataFrame
gdf
Out[15]:
osm_id shop name distance_meters geom
0 4370365690 coffee Kaffeepur 277.208683 POINT (8.51922 47.37178)
1 4515947589 coffee Kaffeewerkstatt 351.997864 POINT (8.5175 47.3707)
2 5297627699 coffee Rent a Barista 797.665201 POINT (8.53152 47.37352)
3 4891858465 coffee Tchibo 1017.406427 POINT (8.535 47.36757)
4 730380415 coffee ViCAFE 1053.386237 POINT (8.53551 47.3676)
5 10967856734 coffee Beanbank 1093.757918 POINT (8.53654 47.36939)
6 10795431343 coffee Bean Bank Coffee & CO 1182.673128 POINT (8.53304 47.37787)
7 11350560294 coffee Nespresso Boutique 1309.342540 POINT (8.53812 47.37478)
In [16]:
# Ensure the GeoDataFrame has the correct projection (EPSG:4326)
if gdf.crs is None:
    gdf.set_crs(epsg=4326, inplace=True)
else:
    pass

# Calculate the mean longitude and latitude for the map center
centroids = gdf.geometry.centroid
lon = centroids.x.mean()
lat = centroids.y.mean()

# Initialize the map
m = folium.Map(location=[lat, lon], 
               zoom_start=16, 
               tiles='ESRIWorldImagery')

# Map settings
folium.GeoJson(
    gdf,
    name='map',
    popup=folium.GeoJsonPopup(fields=['name', 'distance_meters'])
).add_to(m)

folium.LayerControl().add_to(m)

# Plot map
m
Out[16]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Query: Select all roads classified as 'motorway' and create a 5000m buffer around these roads.¶

In [17]:
# Create SQLAlchemy Engine
engine = create_engine(db_connection_url)

# Define SQL query (major roads)
sql = """-- Create buffer around major roads
        SELECT 
            1 as group_id,
            ST_TRANSFORM(ST_UNION(ST_Buffer(p.way::geometry, 5000)), 4326) AS geom
        FROM public.planet_osm_roads AS p
        WHERE
            highway = 'motorway';"""

# Query the database and store the result in a GeoDataFrame
gdf = gpd.GeoDataFrame.from_postgis(sql, engine, geom_col='geom')

# Dispose the engine
engine.dispose()

Show selected features on map¶

In [18]:
# Ensure the GeoDataFrame has the correct projection (EPSG:4326)
if gdf.crs is None:
    gdf.set_crs(epsg=4326, inplace=True)
else:
    pass

# Calculate the mean longitude and latitude for the map center
centroids = gdf.geometry.centroid
lon = centroids.x.mean()
lat = centroids.y.mean()

# Initialize the map
m = folium.Map(location=[lat, lon], 
               zoom_start=9, 
               tiles='EsriWorldTopoMap')

# Map settings
folium.GeoJson(
    gdf,
    name='map'
).add_to(m)

folium.LayerControl().add_to(m)

# Plot map
m
Out[18]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Same ercise but with other type of roads.¶

In [19]:
# Create SQLAlchemy Engine
engine = create_engine(db_connection_url)

# Define SQL query (major roads)
sql = """-- Create buffer around major roads
        SELECT 
            1 as group_id,
            ST_TRANSFORM(ST_UNION(ST_Buffer(p.way::geometry, 2000)), 4326) AS geom
        FROM public.planet_osm_roads AS p
        WHERE
            highway = 'primary';"""

# Query the database and store the result in a GeoDataFrame
primarygdf = gpd.GeoDataFrame.from_postgis(sql, engine, geom_col='geom')

# Dispose the engine
engine.dispose()
In [20]:
# Ensure the GeoDataFrame has the correct projection (EPSG:4326)
if primarygdf.crs is None:
    primarygdf.set_crs(epsg=4326, inplace=True)
else:
    pass

# Calculate the mean longitude and latitude for the map center
centroids = primarygdf.geometry.centroid
lon = centroids.x.mean()
lat = centroids.y.mean()

# Initialize the map
m = folium.Map(location=[lat, lon], 
               zoom_start=9, 
               tiles='EsriWorldTopoMap')



folium.GeoJson(
    primarygdf,
    name='map',
    fill_color='green',
    fillOpacity=0.4,
).add_to(m)


folium.GeoJson(
    gdf,
    name='map',
    fill_color='red',
    fillOpacity=0.4,
).add_to(m)

folium.LayerControl().add_to(m)

# Plot map
m
Out[20]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Jupyter notebook --footer info-- (please always provide this at the end of each notebook)¶

In [21]:
import os
import platform
import socket
from platform import python_version
from datetime import datetime

print('-----------------------------------')
print(os.name.upper())
print(platform.system(), '|', platform.release())
print('Datetime:', datetime.now().strftime("%Y-%m-%d %H:%M:%S"))
print('Python Version:', python_version())
print('-----------------------------------')

-----------------------------------
POSIX
Linux | 6.5.0-1025-azure
Datetime: 2024-10-06 13:34:56
Python Version: 3.12.1
-----------------------------------